Minimal Data Copy for Dense Linear Algebra Factorization

We describe a new result that shows that representing a matrix A as a collection of square blocks will reduce the amount of data reformating required by dense linear algebra factorization algorithms from O(n) to O(n). 1 Description of a Fortran and C Inefficiency for Dense Linear Agebra Factorizations The current Most Commonly Used (MCU) Dense Linear Algebra (DLA) algorithms for serial and SMP processors have a performance inefficiency and hence they give sub-optimal performance. We show that standard Fortran and C two dimensional arrays are the main reason for the inefficiency. We show how to correct these performance inefficiencies by using New Data Structures (NDS) along with so-called kernel routines. The NDS generalizes the current storage layouts for both the Fortran and C programming languages. The BLAS(Basic Linear Algebra Subroutines) were introduced to make the algorithms of DLA performance-portable. However, a relationship exists between the Level 3 BLAS used in most of level 3 factorization routines. This relationship introduces a performance inefficiency in block based factorization algorithms and we will now discuss the Level 3 BLAS, DGEMM (Double precision GEneral Matrix Matrix) to illustrate this fact. In [5, 2] design principles for producing a high performance “Level 3” DGEMM BLAS are given. A key design principle for DGEMM is to partition its matrix operands into submatrices and then call an L1 kernel routine multiple times on its submatrix operands. Another key design principle is to change the data format of the submatrix operands so that each call to the L1 kernel can operate at or near the peak Million FLoating point OPerations per Second (MFlops) rate. This format change and subsequent change back to standard data format is a cause of a performance inefficiency in DGEMM. The DGEMM interface definition requires that its matrix operands be stored as standard Fortran or C two-dimensional arrays. Any DLA factorization algorithm (DLAFA) of a matrix A calls DGEMM multiple times with all its operands being submatrices of A. For each call data copy will be done; the principle inefficiency is therefore multiplied by this number of calls. However, this inefficiency can be eliminated by using the NDS to create a